Code
library(tidyverse)
library(here)
library(tmap)
library(sf)
library(ggExtra)
library(patchwork)library(tidyverse)
library(here)
library(tmap)
library(sf)
library(ggExtra)
library(patchwork)# Load data
ice_area <- read_csv(here('sea_ice_data', 'sibt_areas_v2.csv'))
ice_extent <-read_csv(here('sea_ice_data', 'sibt_extents_v2.csv'))
ice_monthly <- readxl::read_excel(here("sea_ice_data", "Sea_Ice_Index_Monthly_Data_by_Year_G02135_v3.0.xlsx"))
roc_arctic <- readxl::read_excel(here("sea_ice_data", "Sea_Ice_Index_Rates_of_Change_G02135_v3.0.xlsx"))
shapefile <- read_sf(here('ARPA_polygon', 'ARPA_polygon.shp'))
latlong <- readxl::read_xlsx(here('sea_ice_data', 'arctic_regions_latlong.xlsx'))
biomes <- read_csv(here::here('data', 'FireALTEstimatedPairsBurnedUnburned.csv'))
daily <- readxl::read_xlsx(here('sea_ice_data', 'Sea_Ice_Index_Daily_Extent_G02135_v3.0.xlsx'))# Tidy daily ice data
daily_tidy <- daily |>
pivot_longer(cols = c(3:53), names_to = 'year', values_to = 'extent_index') |>
#filter(!is.na(...1)) |>
rename(month = ...1,
day = ...2) |>
mutate(
month = factor(month, levels = month.name, labels = month.abb), # Optional: Convert to abbreviated month names
day = as.numeric(day) # Ensure day is numeric
) |>
fill(month, .direction = "down") |>
mutate(year = as.numeric(str_extract(year, "\\d+"))) |>
filter(year %in% c(2000, 2020))
ggplot(daily_tidy, aes(x = day, y = month, fill = extent_index)) +
geom_tile(color = "white") + # Use white borders between tiles
scale_fill_viridis_c(option = "mako", name = "Extent Index") + # Adjust color scale
facet_wrap(~year) + # Separate heatmaps for each year
theme_minimal() +
labs(
title = "Daily Sea Ice Extent Heatmap Over 20 Year Period",
x = "Day of the Month",
y = "Month",
caption = "Over the 20 year period, sea ice extent is lower on average\n in 2020 compared to 2000 in summer and winter months.\n \nSource: National Oceanic and Atmospheric Administration (NOAA)"
) +
theme(
axis.text.x = element_text(hjust = 1), # Rotate x-axis labels for readability
legend.position = "bottom",
text=element_text(size=15, family="AppleGothic"),
plot.caption = element_text(hjust = .5),)# Rename first column as "YYYYDDD" and the rest based on the first row
colnames(ice_extent) <- ice_extent[1,]
ice_extent <- ice_extent[-1,] # Remove the row used for column names
# Remove the first row (which contains "RegnArea")
ice_area_cleaned <- ice_area[-1, ]
# Convert to long format
extent_tidy_locs <- ice_extent %>%
pivot_longer(cols = 2:18, names_to = "Region", values_to = "ice_extent") %>%
slice(-(1:18)) |> # Fine to start with year 1850
mutate(#YYYDDD = (as.string(YYYYDDD)),
ice_extent = as.numeric(ice_extent)) |> # extent in square kilometers
mutate(year = sub("^(.{4}).*", "\\1", YYYYDDD)) |>
select(-YYYYDDD) |>
group_by(Region, year) |>
summarise(ice_ave = mean(ice_extent)) |>
ungroup() |>
mutate(region = str_trim(Region)) |>
filter(Region != 'Northern_Hemisphere') |>
left_join(latlong)
tidy_monthly_ice <- ice_monthly |>
janitor::clean_names() |>
rename(year = x1) |>
select(-x14) |>
pivot_longer(cols = 2:13,
names_to = "month",
values_to = "ice_extent") |>
mutate(month = str_to_title(month)) |> # Capitalize first letter
mutate(month = factor(month, levels = month.name, labels = month.abb)) ggplot(tidy_monthly_ice, aes(x = month, y = ice_extent, group = year, colour = year)) +
geom_line(size = 1) +
scale_color_viridis_c(option = 'mako') +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = 0.05)) +
theme_classic() +
labs(y = 'Sea Ice Extent Index',
title = 'Seasonal Fluctuations of Sea Ice Extent',
x = ' ',
caption = 'Source: National Oceanic and Atmospheric Administration (NOAA)',
color = 'Year') +
theme(legend.position = 'bottom',
#legend.box.margin = margin(t = 10, r = 10, b = 10, l = 10), # Adjusts legend margins
# plot.margin = margin(t = 15, r = 15, b = 15, l = 15),
legend.key.width = unit(2, "cm"), # Stretches legend keys (adjust as needed)
#legend.spacing.x = unit(0.5, "cm"),
text=element_text(size=15, family="AppleGothic"))biome_plot <- ggplot(biomes) +
geom_area(aes(x = year, y = estDepth, group = distur, fill = distur)) +
scale_fill_manual(values = c('#16425b', '#3a7ca5')) +
# scale_fill_viridis_d(option = 'mako') +
theme_classic() +
# geom_line(data = tidy_monthly_ice, aes(x = year, y =annual), size=2, color = '#023e8a') +
theme_classic() +
labs(x = ' ',
y = 'Estimated Depth (cm)',
title = 'Sea Ice Disturbance and Arctic Permafrost',
subtitle = 'Increasing Active Layer Thickness (ALT) in the Arctic Tundra\nmelts the permafrost layer more sporadically due \nto rising global temperatures',
caption = 'ALT refers to the thickness of the layer above\npermafrost that freezes and thaws seasonally. \nData Source: Arctic Data Center',
fill = 'Disturbance') +
theme(text=element_text(size=13, family="AppleGothic"),
# legend.title = element_text(size = 12, hjust = .5), #sets legend title size
legend.position = c(.25, .25), #sets legend to the top
# legend.position = 'bottom',
legend.text = element_text(size = 10),
plot.caption = element_text(hjust = .5),
panel.border = element_rect(color = "black", fill = NA, size = .5)) +
scale_x_continuous(expand = c(0,0), breaks = scales::pretty_breaks(), position = 'top') + # Move x-axis to the top
scale_y_reverse(expand = c(0,0))Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
#biome_plot
biomes$permaExtent <- factor(biomes$permaExtent,
levels = c('C', 'D', 'S'),
labels = c('Continuous', 'Discontinuous', 'Sporadic'))
p2 <- ggplot(biomes) +
geom_density(aes(y = estDepth, fill = permaExtent), alpha = 0.5) + # Swap x → y
scale_y_reverse() + # Flip so 0 is at the top
scale_fill_viridis_d(option = 'mako') +
theme_void() + # Remove extra elements for a clean look
labs(fill = 'Permafrost Extent') +
theme(text=element_text(size=12, family="AppleGothic"))
#legend.position = 'bottom')
#p2
# Combine the two plots
combined_viz <- biome_plot + p2 + plot_layout(widths = c(4, 1)) # Main plot wider than density plot
combined_viz#ggsave(plot = combined_viz, filename = 'arctic-permafrost.png', width = 8, height = 6)I still plan on pursuing option 1, the infographic.
When do we see the greatest reduction or increase in sea ice extent? How has sea ice changed over many years? How does decreasing sea ice extent impact permafrost layers?
My first two questions are largely the same to my initial proposal, but I decided to ask the last question about permafrost because I think it is an interesting component to supplement these trends! I found a cool dataset from the Arctic Data Center that I incorporated.
I will be using year, month, date, average layer thickness depth, and other categorical permafrost variables to help answer my questions. Much of the original NOAA data did not have many categorical variables like permafrost extent and disturbance, so this dataset from the Arctic Data Center have some really interesting observations that classify their measurements with categories. Since there is seasonal fluctuations, I want to show the daily ice extent along with the monthly ice extent curve, portrayed in different ways ,a heatmap and a line plot. My permafrost plot is an area plot that looks like icicles, and I want to show density curves of different permafrost extent indicators based on the ALT depth. This will answer my last question, and the other temporal plots will answer the first two questions.
Sample Data Viz 1: Derek Watkins - ‘Yearly fluctuations in are of Arctic covered by ice’ The sample data viz above was shown in class and it is so cool! I love that they highlighted important lines and pointed out specific points. I think we are using the same data so it is definitely possible, but I like the annotation elements to the text. I am curious if there is a way to annotate things to frame it with permafrost. For example, as ice extent is on the decline, perhaps there is an event that I can note that is not already on the sample visualization. I think it is now a question of which lines are the most important to highlight.
Sample Data Viz #2: Philippe Rekacewicz - Where is permafrost located?
I really want to add a geospatial component to my infographic. I have a shapefile of the arctic circle, but it is just a single polygon and I definitely need a refresher on geospatial analyses. I like how they highlgihted the permafrost layers through a map, because my data uses the same categories. I also have another data set with locations of sea ice area, and I manually went into google earth and extracted coordinates, so I am curious if I am able to make something like this combining different sites, or perhaps just highlighting certain points in the arctic circle with a bubble chart where the size of the bubble corresponds with lost ice area.